% March, 2012
%
% This program sets up a simple external scattering problem for Helmholtz.

function res = getnpts(nobj,kh)

global anim;
global nobj0;
global level;
global xctot;
global yctot;
global Rtot;
global isave;

anim = false;
check_sol = true;

% number of points per wavelength
ppl = 10;

isave = 0;

nobj0 = nobj;

% Angle of the incoming wave

theta_in  = 0;

n = log(nobj) / log(2);

n1 = 2^(floor(n/2));
n2 = nobj / n1;

centers = [];
for i2=1:n2
  for i1=1:n1
    centers = [centers;[i1*4 + mod(i2,2)*2,i2*2.5]];
  end
end
nobjsx = n1;
nobjsy = n2;

n   = 100;  % The number of points on each contour.
acc = 1e-10;

fprintf(1,'\n=== Running a test on %d contours.\n',nobj)
fprintf(1,'    Frequency : %g \n\n',kh)

% Set the geometry.
contours = cell(nobj,1);
matsA    = cell(nobj,1);

R        = zeros(nobj,1);
xc       = zeros(nobj,1);
yc       = zeros(nobj,1);
ntot = n * nobj;

hmax_phy = 0;
figure(1); hold on
Ctot = [];
for i=1:nobj
  angle       = 2*pi*rand();
  k           = floor(6*rand());
  [C,len,xin,tmp] = get_geometry_star(n,centers(i,:).',k,angle);
  h_phy       = max(sqrt(C(2,:).^2 + C(5,:).^2));    
  hmax_phy    = max(hmax_phy,h_phy);
  contours{i} = C;
  xc(i)       = mean(C(1,:));
  yc(i)       = mean(C(4,:));
  R(i)        = max(sqrt((C(1,:) - xc(i)).^2 + (C(4,:) - yc(i)).^2 ));
  plot(C(1,:) ,C(4,:) ,'k')
  Ctot        = [Ctot,C];          
end
xctot       = mean(Ctot(1,:));
yctot       = mean(Ctot(4,:));
Rtot        = max(sqrt((Ctot(1,:) - xctot).^2 + (Ctot(4,:) - yctot).^2 ...
                       ));

lambda = 2*pi/kh;
% plot([xctot - .5*lambda*cos(theta_in),xctot + .5*lambda*cos(theta_in)],...
%      [yctot - Rtot - .5*lambda*sin(theta_in),yctot - Rtot + .5*lambda*sin(theta_in)],...
%      'r','linew',2)



% Check the size of the discretization vs the wavelength
fprintf(1,'h / lambda = %g\n', kh*hmax_phy/(2*pi))

if kh*hmax_phy/(2*pi) > 1./10
  error('mesh too coarse')
end

fprintf(1,'R / lambda = %g\n', kh*R(1)/(2*pi))

% Incoming wave
% Incoming potential

if (check_sol)
  theta_zzz = linspace(0,2*pi,500).';
  xxx       = xctot + 2.5*Rtot*cos(theta_zzz);
  yyy       = yctot + 2.5*Rtot*sin(theta_zzz);
  zzz = [reshape(xxx,1,numel(xxx));...
	 reshape(yyy,1,numel(yyy))];
  figure(1)
  plot(zzz(1,:),zzz(2,:),'r')
  matA   = zeros(ntot);
end

alabel = gca();
set(alabel,'FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
axis off
axis equal
print(gcf, '-dpng', 'config.png');

for i=1:nobj
  i0 = (i-1)*n + 1; i1 = i0 + n - 1;
  matAii = get_A_single_diag(contours{i},1:n,kh);
  matsA{i}          = matAii;
  if (check_sol)
    for j=1:nobj
      i0 = (i-1)*n + 1; i1 = i0 + n - 1;
      j0 = (j-1)*n + 1; j1 = j0 + n - 1;
      if (i == j)
	matA(i0:i1,j0:j1) = matAii; 
      else
	matA(i0:i1,j0:j1) = get_A_offd_noquad(Ctot,i0:i1,j0:j1,kh);
      end
    end
  end
end

if (check_sol)
  v     = -incoming_wave(kh,theta_in,Ctot(1,:),Ctot(4,:)); v = v.';
  sigma = matA\v;
  phi_scattered_ref  = evalpot(zzz,Ctot,sigma,kh);
end


level = 1;
res = [[1;ntot]];
while nobj>=2
  fprintf(1,'\n\n Level %d \n',level)
  level = level + 1;
  
  [contours_up,matsA_up,nobjsx,nobjsy] = factor_level(nobjsx,nobjsy,contours, ...
						  matsA,kh,acc,ppl);        
  contours = contours_up;
  matsA    = matsA_up;
  nobj     = nobjsx*nobjsy;
  
  npts = 0;
  for i=1:nobj
    Ctmp = contours{i};
    npts = npts + size(Ctmp,2);
  end
  res = [res,[level;npts]];
end

[D,h_D]          = get_circle(contours,1,kh,ppl,1.);
[Is,Cs,V,matP,D] = get_skeleton(kh,contours{1},matsA{1},acc,D,h_D);


npts =size(Cs,2);
res = [res,[level+1;npts]];

if (check_sol)
  v     = -incoming_wave(kh,theta_in,Cs(1,:),Cs(4,:)); v = v.';
  % Induced charge distribution
  sigma = matP*v;
  % Potential created on D
  phi_scattered  = evalpot(zzz,Cs,sigma,kh);
  err = max(abs(phi_scattered - phi_scattered_ref));
  fprintf(1,'\n\n ERROR :  %g \n',err)
end


figure(); hold all
plot(res(1,:)-1,res(2,:),'ko-','markersize',20)
xlabel('level','FontName','cmr10','FontSize',15)
ylabel('number of discretization points','FontName','cmr10','FontSize',15)
alabel = gca();
set(alabel,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
print(gcf, '-dpng', 'npts.png');

return
